 function x = read_ncat(varargin)
%function x = read_ncat(options)
%
% read in a slice of the NCAT phantom
% and assign it attenuation coefficients in inverse mm units.
%
% options
%	'nx'		desired size
%	'ny'
%	'mu'	[5,1]	desired image intensity values
%			0 background, 1 lung, 2 body, 3 spine, 4 ribs
%	'ddir'		data directory
% out
%	x	[nx,ny]	ncat image
%

if ~nargout, help(mfilename), read_ncat_test, error(mfilename), end

arg.nx = 1024;
arg.ny = [];
arg.mu = [0:4];
arg.ddir = '';
arg.file = 'ncat,1024,slice500,0to4.fld';
arg = vararg_pair(arg, varargin);
if isempty(arg.ny), arg.ny = arg.nx; end

x = read_ncat_do(arg.ddir, arg.file, arg.nx, arg.ny, arg.mu);


%
% read_ncat_do()
%
function x = read_ncat_do(ddir, file, nx, ny, mu)

% guess .../data directory by looking parallel to '.../transmission' directory
if ~isvar('ddir') | isempty(ddir)
	t = path_find_dir([filesep 'transmission']);
	ddir = strrep(t, 'transmission', 'data');
end

if ~exist(ddir, 'dir')
	warning(sprintf('cannot find data directory %s', ddir))
	error(sprintf('edit path in %s.m for your installation', mfilename))
end

file = [ddir filesep 'ncat,1024,slice500,0to4.fld'];
x = fld_read(file);

x = mu(1+x); % map indices to desired image intensity values

x = phantom_resize(x, nx, ny);


function read_ncat_test
x = read_ncat('nx', 512, 'ny', 450, 'mu', [0 0.05 0.2 0.3 0.4]/0.2*1000);
im(x, 'ncat phantom'), cbar
